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Abstract 

In  this  paper  we  describe  a  new  3rd-order  algorithm  for  solving  the  transport  equa¬ 
tions  of  plasma  in  highly  non-equilibrium  conditions.  The  plasma  is  described  as  a 
two-temperature,  single  fluid  with  the  kinetics  of  collisional  and  radiative  excitation 
and  ionization,  and  reverse  processes.  This  Collisional-Radiative  model  is  currently 
limited  to  atomic  plasma  and  does  not  include  radiative  transport.  We  describe  in 
detail  some  special  techniques  for  level  grouping,  scale  separation  of  slow  (trans¬ 
ported)  and  fast  (quasi-steady-state)  level  kinetics,  and  a  non-linear  transformation 
of  the  transported  equations  of  the  electronic  levels  to  achieve  the  desired  accuracy. 
The  implementation  and  testing  of  the  various  coupling  and  relaxation  processes  are 
described.  The  fluid  transport  is  computed  using  a  3rd-order  variant  of  the  MP5 
monotonicity-preserving  upwind  advection  scheme.  The  code  is  implemented  in  Java 
and  parallelized  through  domain  decomposition  and  hierarchical  multi-threading;  ap¬ 
proach  and  performance  are  also  briefly  discussed.  The  numerical  model  is  validated 
on  various  standard  test  cases,  and  applied  to  the  simulation  of  ionizing  shock  front 
propagation  in  Argon.  This  problem  shows  a  high  sensitivity  to  the  kinetics  ladder 
of  ionization  and  population  of  the  excited  states,  leading  to  fluctuations  of  the  lo¬ 
cation  of  the  electron  avalanche  at  the  end  of  the  induction  zone  behind  the  shock. 
We  show  that  the  collisional-radiative  kinetics  can  reproduce  the  corrugations  of  the 
shock  front  observed  in  the  experiment. 


Distribution  A:  Approved  for  public  release;  distribution  unlimited. 
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Nomenclature 


frozen  sound  speed 
total  energy  density 
electronic  thermal  energy  density 
heavy  particle  thermal  energy  density 
specific  stagnation  enthalpy 
shock  Mach  number 
number  density 
total  thermal  pressure 
electronic  thermal  pressure 
heavy  particle  thermal  pressure 
electronic  temperature 
heavy  particle  temperature 
mean  thermal  speed 
heavy  particle  temperature 
electron  temperature 
momentum  density 
mass-average  velocity 
p  mass  density 

a  ionization  fraction 

£  relaxation  length 

a  collision  cross  section 
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1  Introduction 

First  introduced  by  Bates  et.  al.  [2],  Collisional-Radiative  (cr)  models  have  been  exten¬ 
sively  used  in  calculating  population  distribution  of  energy  states  as  well  as  obtaining 
information  about  the  relative  significance  of  the  various  physical  processes.  Here  we  in¬ 
troduce  a  CR  model  coupled  with  a  high-order  convective  scheme  to  model  nonequilibrium 
ionization  kinetics  in  high  speed  gas  flows.  Only  the  first  few  excited  levels  are  included 
in  our  model  which  are  convected  as  individual  fluid  species.  It  is  expected  that  inclusion 
of  only  the  first  few  excited  levels  are  sufficient  to  capture  and  resolve  the  shock  struc¬ 
ture.  With  this  formulation  we  are  potentially  able  to  determine  the  validity  of  these 
assumptions  in  various  plasma  regimes,  and  estimate  their  significance.  This  particular 
modeling  would  be  useful  when  examining  the  transition  of  a  plasma  from  a  collisional  to 
a  radiation-dominated  regime,  or  when  looking  in  detail  at  the  shock  layers  and  relaxation 
regions.  For  the  application  presented  here,  we  assume  the  following: 

•  the  plasma  is  neutral,  and  all  components  have  a  single  velocity  (single  fluid  approx¬ 
imation). 

•  all  components  are  well  represented  by  a  Maxwell-Boltzmann  distribution  of  the 
velocity. 

•  electromagnetic  effects  are  not  considered. 

In  the  present  work,  the  propagation  of  a  strong  ionizing  shock  wave  in  Argon  is  studied, 
and  used  as  a  first  step  in  the  validation  of  our  code. 


2  Transport 

We  consider  a  neutral,  single-fluid,  two-temperature,  nonequilibrium  plasma.  In  a  finite- 
volume  framework,  the  hyperbolic  system  of  equations  has  the  form: 

^  f  Q  dV  +  n-F  dS  =  [  ndV  (1) 

dt  Jv{t)  Js(t )  Jv{t) 
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being  the  total  energy  density  of  the  plasma  and  HQ  =  E  +  p  the  stagnation  enthalpy 
density.  In  order  to  maintain  the  hyperbolic  property  of  all  convected  terms,  the  electron 
energy  equation  is  formulated  using  the  electronic  entropy  function  se  =  pe/ Ple  (cf.  [3]). 
The  inclusion  of  p  as  opposed  to  pe  ensures  that  the  function  is  well-defined  for  all  ionization 
fractions.  This  leads  to  the  following  differential  relations  between  the  conservative  and 
primitive  variables: 
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£sdps - - 
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for  which  we  have  used  the  ideal  gas  equation  of  state  p  =  ^  psRsTs  as  well  as 


lh  -  1 


X)  PsRs 

s^e 

Ps^v,s 

s^e 


(6) 


Notice  that  the  electron  component  contributes  to  the  total  plasma  pressure,  p  = 
Ph  +  pe-  The  electron  energy  equation  is  linear  and  is  responsible  for  the  heating  of 
electrons  through  adiabatic  compression.  This  final  form  of  the  electron  energy  convection 
can  be  derived  from  the  reduction  of  the  two-fluid  system  of  equations  to  the  single  fluid 
system. 

The  term  Cl  contains  the  source  terms  for  the  CR  model.  The  Qa,p  terms  represent  the 
energy  exchange  terms  between  the  electrons  and  various  plasma  components,  a,  (3.  These 
terms  are  anti-symmetric  in  a,  /?  (due  to  energy  conservation)  and  their  sum  vanishes. 
The  to  are  source/sink  terms  for  specie  densities  due  to  ’chemical  reactions’  and  radiative 
processes,  and  Q  are  the  corresponding  change  in  energy  for  a  given  plasma  component 
due  to  ’chemical  reactions’  and  radiative  processes.  The  term  ’chemical  reaction’  refers 
to  any  collisional  process  that  changes  the  number  or  type  of  species  involved.  We  can 
separate  between  the  collisional  and  radiative  parts,  ClfR  =  Clf  +  ClR,  and  since  energy 
conservation  requires  that  ]P.S  Clf  =  0,  the  final  sum  on  the  R.H.S.  is  composed  of  the 
radiative  terms  only. 

To  this  system  of  equations,  one  should  add  fluxes  due  to  viscous  and  dissipative 
processes,  in  particular  electron  heat  conduction,  which  can  be  a  very  rapid  process.  Our 
scheme  also  solves  for  this  process. 

The  above  system  is  solved  using  the  operator-splitting  approach  under  the  Quasi- 
Steady-State  (qss)  approximation.  Using  this  approach,  the  convective  terms  are  solved 
independently  of  the  source  terms  via 


dt  +v-^n 


,A»  =  0 


s 


(7) 
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which  is  then  explicitly  integrated  in  time  using  the  second-order  Adams-Bashforth  scheme. 
The  source  terms,  on  the  other  hand,  are  treated  implicitly,  requiring  a  computation  of  a 
dense  Jacobian  matrix.  Upwind  fluxes  at  the  cell  interfaces  are  evaluated  using  the  HLLE 
Riemann  solver  [4].  High-order  spatial  resolution  is  achieved  via  parabolic  interpolation 
of  the  left  and  right  states  at  the  cell  interfaces  via 


Ql  =  +  5 qj  -  qj+1) 

QR  =  ^(2gj+i  +  5 qj  -  qj- i), 


(8) 

(9) 


where  q  are  the  conserved  quantities.  The  above  interpolation  is  third-order  accurate  [8]. 
Strong  nonlinear  waves  require  limiting  which  modifies  the  left  and  right  states  according 
to 


qL  <-  median(gL,  qj ,  qMP)  (10) 

with  qAIP  =  gj+minmod^+i— qj,  ot(qj—  being  the  monotonicity-preserving  limit  [7]. 

The  value  of  a  is  taken  to  be  4.  The  right  state  is  found  from  symmetry. 

3  Collisional-radiative  model 

We  have  considered  the  following  processes  in  our  CR  model, 
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Fji 

(11) 

Ar{i)  +  Ar(  1)  N  N  Ar{j )  +  Ar(  1) 

Lji 

(12) 

Ar(i )  +  e  Ar+  +  e  +  e 

Oi 

(13) 

Ar{i)  +  Ar(  1)  Ar+  +  e  +  Ar(  1) 

Wi 

(14) 
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where  the  rate  coefficients  defined  in  [12]  are  summarized  in  table  2. 

In  our  calculations  we  have  considered  only  the  excited  levels  of  the  3p54s  manifold 
(cf.  table  3).  This  implies  that  ionization  and  recombination  should  proceed  from  and  to 
only  these  low-lying  levels;  although  levels  beyond  this  manifold  are  more  that  leU  away 
from  the  ionization  limit,  the  combination  of  small  energy  gap  and  large  cross-section 
makes  the  ionization  from  these  levels  extremely  rapid,  certainly  with  time  scales  much 
below  the  time  resolution  needed  in  our  computations. 
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Table  2:  Rate  coefficients  for  collisional-radiative  model. 


Rate 

Coefficient 

Process 

Cnm 

collisional  excitation  by  electrons 

Knm 

collisional  excitation  by  ground  state  atoms 

Fmn 

collisional  de-excitation  by  electrons 

Lmn 

collisional  de-excitation  by  ground  state  atoms 

Sm  5 

collisional  ionization 

Om,  Wm 

three-body  recombination 

Rm 

radiative  recombination 

A-mn 

transition  probability  (Einstein  coefficient) 

A-mn 

bound-bound  optical  escape  factor 

A  771 

bound-free  optical  escape  factor 

Bm 

three-body  collision  of  metastable  and  ground  state  atoms 

Table  3:  Levels  considered  in  collisional-radiative  model. 
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Table  4:  Parameters  related  to  atom-atom  processes. 
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3.1  Excitation  and  ionization  by  heavy  particle  collisions 

Since  we  are  considering  a  shock  propagating  into  a  non-ionized  gas,  the  inelastic  atom- 
atom  processes  are  of  critical  importance.  Atom-atom  excitation  is  the  most  important 
controlling  process  for  preliminary  ionization  as  well  as  the  overall  relaxation  length.  Un¬ 
fortunately,  very  little  experimental  data  is  available  for  the  cross  sections  of  such  pro¬ 
cesses.  However,  due  to  the  relatively  low  energy  of  the  flow  under  consideration,  linear 
approximations  are  sufficient.  The  cross  sections  for  excitation  from  ground  state 


=Pij(£-£ij)  (17) 

as  well  as  for  the  inner  3p54s  manifold  transitions 

<4(e)  =  P!j (g  2.26  ^  (18) 

are  taken  from  [1]  with  the  parameter  /3*?-  given  in  table  4.  For  direct  ionization  from 
ground  state  Argon  we  have 

al-Ar(e)  =  l-8  X  1(U25(£  -  Ej)1'3  [m2].  (19) 

The  corresponding  rate  equations  are 
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Table  5:  Parameters  for  the  calculation  of  <r®-.  A,  S ,  and  P  represent  allowed,  spin,  and 
parity  forbidden  transitions  respectively. 
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3.2  Excitation  and  ionization  by  electron  impact 

Excitation  rates  for  electron  impact  processes  are  given  by 


(8kBTe\1/2  [°° 
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where  the  corresponding  cross  sections  are  given  by 
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The  necessary  parameters,  /mna^n,  or  a^n,  are  given  in  table  5. 
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3.3  Photorecombination 


In  the  absence  of  a  third  body,  the  energy  is  liberated  as  radiation.  Photorecombination 
is  a  significant  loss  mechanism  and  plays  an  important  role  in  radiative  cooling.  The  rate 
for  this  process  is  given  by 
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with  cross  sections 
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3.4  Bremsstrahlung  emission 

Radiative  cooling  is  taken  into  account  via  bound-free  Bremsstrahlung  radiation.  This  is 
given  by  the  Kramers-Unsold  formula: 
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327 r  /  27t 

3  \3  mekBTe 
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(31) 


3.5  Elastic  collisions 

Elastic  collisions  are  incorporated  into  the  CR  implicit  solver  as  well.  This  strong  coupling 
permits  more  accurate  and  stable  calculations.  The  energy  transfer  between  electrons  and 
heavy  particles  is  computed  as 

9Ee  _  3^  ^  2 me  frr  rr  /0r,\ 

oi  —  V^h  TejCM.e 

( 'JL  jL  777/ 

where  the  collision  rate  is  given  by 


For  electron-Argon  collisions,  we  have  used  the  theoretical  cross-sections  given  in  ta¬ 
ble  3.5  as  computed  by  McEachran  and  Stauffer  [6].  The  corresponding  electron-Argon 
and  electron-Argon+  collision  rates  are  given  in  figures  3.5  and  3.5  respectively. 


4  Results  and  discussion 

The  code  was  applied  to  the  case  of  a  strong  ionizing  shock  in  Argon,  in  an  attempt  to 
reproduce  a  series  of  experiments  by  Glass  &  Liu  [5] .  The  experimental  cases  chosen  were 
for  a  shock  Mach  number  of  14.7.  The  typical  shock  structure  contains  some  important 
observable  quantities  which  are  defined  here: 

•  Ms  is  the  instantaneous  shock  Mach  number. 

•  G  is  the  relaxation  length,  i.e.  the  distance  from  the  shock  front  to  the  location  of 
the  peak  ionization. 


Table  6:  Cross  sections  for  e-Ar  elastic  collisions. 


E  [eV] 

cr  x  1020  [m2] 

0.01 

4.4679 

0.03 

2.9180 

0.05 

2.1193 

0.07 

1.6052 

0.09 

1.2481 

0.13 

0.7972 

0.17 

0.5428 

0.19 

0.4621 

0.21 

0.4002 

0.23 

0.3548 

0.25 

0.3242 

0.29 

0.2934 

0.31 

0.2898 

0.32 

0.2904 

0.41 

0.3504 

0.51 

0.4756 

0.61 

0.6403 

0.71 

0.8240 

0.81 

1.0184 

0.91 

1.2176 

1.00 

1.4002 

1.50 

2.4307 

2.00 

3.4680 

3.00 

5.5581 

4.00 

7.7317 

5.00 

10.0665 

7.50 

16.7176 

10.00 

22.4036 

10 
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Figure  1:  Electron-Ar  elastic  collision  integral. 
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Figure  2:  Electron-Ar+  Coulomb  collision  integral. 
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Figure  3:  Results  of  collisional-radiative  kinetics  only-electron  and  heavy  particle  temper¬ 
atures  as  well  as  ionization  fraction. 

•  a*  is  the  peak  ionization  fraction. 

•  r*  shock  fluctuation  periodicity. 

Results  of  the  stand-alone  collisional-radiative  kinetics  solver  are  shown  in  figures  3 
and  4.  The  simulated  plasma  conditions  are  those  corresponding  to  a  shock  propagating 
at  Mach  14.7  into  an  argon  gas  with  Th  =  300  K  and  p  =  544  Pa.  Evident  is  the  initial 
plasma  formation  and  thermalization  between  electrons  and  heavy  particles  followed  by 
the  electron  avalanche  and  radiative  cooling  of  the  gas.  The  electron  temperature  rises 
rapidly  to  a  plateau,  then  rises  very  slowly,  until  a  maximum  which  corresponds  to  the 
avalanche  position,  then  decreases  until  equilibrium  with  the  heavy  particles  is  achieved. 
The  resulting  time-dependent  profile  is  indicative  of  a  steady-state  shock  structure  under 
the  same  conditions. 

For  coupling  between  the  CR  model  and  the  fluid  transport  scheme,  a  shock  was  im¬ 
pulsively  started  by  modeling  a  low-pressure  Argon  gas  traveling  at  high  velocity  from 
left  to  right  and  reflecting  off  a  wall  at  the  edge  of  our  grid  system.  The  instantaneous 
measurements  of  the  shock  Mach  number  were  accordingly  transformed  to  the  rest  frame 
of  the  gas.  Since  the  calculation  is  unsteady,  we  were  able  to  observe  the  initial  formation 
of  the  plasma,  and  its  coupling  to  the  shock  front;  prior  to  ionization,  all  of  the  energy  in 
invested  in  the  translational  modes,  and  the  shock  speed  is  much  greater  than  the  final, 
stable  value,  as  will  be  shown  shortly. 
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Figure  4:  Results  of  collisional-radiative  kinetics  only-specie  number  densities. 
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Figure  5:  Temperature  profiles  at  three  different  times  as  shock  propagates  away  from  the 
wall.  Snapshots  from  right  to  left  are  taken  at  8.4,  35,  and  82  /isec  with  corresponding 
Mach  numbers  of  approximately  16,  14,  and  15,  respectively. 
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Figure  6:  Shock  Mach  number  as  a  function  of  distance. 


Figure  5  shows  the  shock  structure  obtained  for  the  case  of  a  shock  traveling  with 
an  initial  Mach  number  of  16  into  Argon  with  Th  =  300  K  and  p  =  544  Pa,  after 
3  distinct  distances  from  the  reflecting  wall.  These  oscillations  of  the  peak  ionization, 
relaxation  length,  and  shock  Mach  number  were  identified  as  true  physical  phenomena.  An 
explanation  for  the  coupling  mechanism  can  be  devised  as  follows:  a  pressure  fluctuation 
travels  towards  the  shock,  and  changes  its  velocity.  As  the  shock  strength  is  affected  by 
the  pressure  fluctuation,  the  post-shock  condition  is  slightly  altered.  A  fluid  element  will 
go  through  this  slightly  different  environment  until  the  ionization  avalanche.  The  overall 
periodicity  should  therefore  be: 


T* 


4 


1 


a2  -  u2 


(34) 


For  a2  =  2890  m/s  and  u2  =  3841  m/s,  this  gives  a  value  of  24  psec.  Roughly  one  period 
of  the  oscillatory  behavior  of  the  shock  Mach  number  and  relaxation  length  is  shown  in 
figure  6.  The  initial  shock  Mach  number  of  16  decays  to  below  14  before  rising  again  to 
about  15. 
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5  Concluding  remarks 

A  collisional  radiative  database  has  been  compiled  and  successfully  utilized  in  the  simula¬ 
tion  of  non-equilibrium  ionizing  flows  in  Argon.  The  preliminary  conclusions  of  this  study 
are  as  follows: 

•  The  computed  rates  have  been  successfully  applied  to  flows  from  300  —  25,  OOOif. 

•  Inclusion  of  ground  state  argon  with  the  first  four  levels  is  sufficient  to  accurately 
reproduce  the  correct  shock  structure  and  relaxation  length. 

•  Fluctuations  in  the  shock  velocity,  coupled  to  the  fluctuations  in  the  electron  den¬ 
sity  at  the  avalanche,  have  been  identified  and  investigated.  These  fluctuations  are 
physical  and  have  a  definite  periodic  character.  A  coupling  mechanism  through  the 
propagation  of  pressure  waves  in  the  relaxation  region  has  been  suggested. 

The  time-accurate  code  used  here  clearly  is  a  powerful  tool  that  can  be  used  to  eventu¬ 
ally  supplement  the  experimental  diagnostic  techniques  and  provide  a  test-bed  for  physical 
models.  The  computational  cost  is  essentially  dominated  by  the  nonequilibrium  kinetics. 
For  more  complex  systems,  notably  shock  layers  around  re-entry  vehicles,  the  appropriate 
method  would  consist  of  successive  approximations  to  the  steady  state  with  sub-cycling. 
The  greatest  advantage  of  the  present  capability  is  the  ability  to  study  unsteady  phenom¬ 
ena,  and  this  will  be  emphasized  in  the  future. 


6  Appendix 

For  completeness,  we  give  the  eigensystem  for  the  two-temperature  model. 
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where 

U±  =  Yj  ysEPa  +  —EPe  +  a2 Ep  ±  -EV6 

'  p  P 

and  the  frozen  sound  speed: 

a2  =  ysPP*  +  [{E+p-pu-  u  )PE  +  SePs  J  /p 


(35) 


(36) 
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